##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)

##########################################################################################

option_list <- list(
    make_option(c("--gene"), type = "character") ,
    make_option(c("--mut_rate_gene_file"), type = "character") ,
    make_option(c("--mut_rate_point_file"), type = "character") ,
    make_option(c("--images_path"), type = "character")
)

if(1!=1){
    
    gene <- "TP53"
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    mut_rate_gene_file <- paste(work_dir,"/images/mutRate/MutRate.tsv",sep="")
    mut_rate_point_file <- paste(work_dir,"/images/mutRate/MutRate.RecurrentPoint.tsv",sep="")
	images_path <- paste(work_dir,"/images/mutRatePlot",sep="")

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gene <- opt$gene
mut_rate_gene_file <- opt$mut_rate_gene_file
mut_rate_point_file <- opt$mut_rate_point_file
images_path <- opt$images_path

dir.create(images_path , recursive = T)

###########################################################################################

col <- c(
  brewer.pal(9,"YlGnBu")[6],
  rgb(234,106,79,alpha=255,maxColorValue=255),
  rgb(203,24,30,alpha=255,maxColorValue=255),
  rgb(255,0,0,alpha=255,maxColorValue=255)
  )

names(col) <- c("IM" , "IGC" , "DGC" , "GC")

###########################################################################################

dat_mutRateGene <- data.frame(fread( mut_rate_gene_file ))
dat_mutRatePoint <- data.frame(fread( mut_rate_point_file ))

###########################################################################################

dat_mutRateGene <- subset(dat_mutRateGene , Hugo_Symbol==gene)
dat_mutRatePoint <- subset(dat_mutRatePoint , Hugo_Symbol==gene)
dat_mutRatePoint <- subset( dat_mutRatePoint , vid %in% unique(subset(dat_mutRatePoint , MutNum > 5)$vid) )
dat_mutRatePoint$Hugo_Symbol <- dat_mutRatePoint$vid
dat_mutRatePoint <- dat_mutRatePoint[,-3]

dat_plot <- rbind( dat_mutRateGene , dat_mutRatePoint )
dat_plot$Class <- factor( dat_plot$Class , levels = c("IM" , "GC" , "IGC" , "DGC") , order = T )
dat_plot$value_text <- paste0( round(dat_plot$MutRate , 2) * 100 , "%") 

###########################################################################################

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("p == ",down," %*% 10","^",up)
    return(text)
}


###########################################################################################
## 计算P值

dat_plot$p.value = ""
dat_plot$p_text = ""

result <- c()
for(geneN in unique(dat_plot$Hugo_Symbol)){
    for(from in unique(dat_plot$From)){

        print(geneN)

        tmp_1 <- subset( dat_plot , Hugo_Symbol == geneN & Class %in% c("GC") & From == from )
        tmp_2 <- subset( dat_plot , Hugo_Symbol == geneN & Class %in% c("IM") & From == from )

        if(nrow(tmp_1) != 0 | nrow(tmp_2) != 0){
            if(nrow(tmp_1)==0){
                tmp_1 <- tmp_2
                tmp_1$MutNum <- 0
                tmp_1$MutRate <- 0
                tmp_1$value_text <- ""
                tmp_1$Class <- "GC"
            }

            if(nrow(tmp_2)==0){
                tmp_2 <- tmp_1
                tmp_2$MutNum <- 0
                tmp_2$MutRate <- 0
                tmp_2$value_text <- ""
                tmp_2$Class <- "IM"
            }


            tmp <- rbind( tmp_1 , tmp_2 )

            tmp_fisher <- matrix(c(tmp$MutNum , tmp$SampleNum - tmp$MutNum) , ncol = 2)
            p <- fisher.test(tmp_fisher)$p.value

            if( p < 0.001 ){
                p_text <- trans(p)
            }else{
                p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
            }

            tmp$p_text <- p_text

            result <- rbind( result , tmp )
        }
    }
}

for(geneN in unique(dat_plot$Hugo_Symbol)){
    for(from in unique(dat_plot$From)){

        print(geneN)

        tmp_1 <- subset( dat_plot , Hugo_Symbol == geneN & Class %in% c("IGC") & From == from )
        tmp_2 <- subset( dat_plot , Hugo_Symbol == geneN & Class %in% c("DGC") & From == from )

        if(nrow(tmp_1) != 0 | nrow(tmp_2) != 0){
            if(nrow(tmp_1)==0){
                tmp_1 <- tmp_2
                tmp_1$MutNum <- 0
                tmp_1$MutRate <- 0
                tmp_1$value_text <- ""
                tmp_1$Class <- "IGC"
            }

            if(nrow(tmp_2)==0){
                tmp_2 <- tmp_1
                tmp_2$MutNum <- 0
                tmp_2$MutRate <- 0
                tmp_2$value_text <- ""
                tmp_2$Class <- "DGC"
            }

            tmp <- rbind( tmp_1 , tmp_2 )

            tmp_fisher <- matrix(c(tmp$MutNum , tmp$SampleNum - tmp$MutNum) , ncol = 2)
            p <- fisher.test(tmp_fisher)$p.value

            if( p < 0.001 ){
                p_text <- trans(p)
            }else{
                p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
            }

            tmp$p_text <- p_text
            result <- rbind( result , tmp )
        }
    }
}

dat_plot <- result

dat_plot$Hugo_Symbol <- factor( dat_plot$Hugo_Symbol , 
    levels = unique(dat_plot[order(dat_plot$MutRate , decreasing=T),"Hugo_Symbol"]) , order = T)

###########################################################################################

for( gene in unique(dat_plot$Hugo_Symbol) ){
    dat_plot2 <- subset( dat_plot , Hugo_Symbol == gene & From == "All" )
    dat_plot2$Type <- ""
    dat_plot2$Type[dat_plot2$Class %in% c("IM" , "GC")] <- "IM vs GC"
    dat_plot2$Type[dat_plot2$Class %in% c("IGC" , "DGC")] <- "IGC vs DGC"
    dat_plot2$Type <- factor( dat_plot2$Type , levels = c("IM vs GC" , "IGC vs DGC") )
    dat_plot2$Class_num <- paste0( dat_plot2$Class , "(" , dat_plot2$SampleNum , ")" )  
    dat_plot2$Class_num <- factor( dat_plot2$Class_num , levels = unique(dat_plot2$Class_num)[order(unique(dat_plot2$Class_num) , decreasing=T)] )  

    ###########################################################################################
    ## IGC vs DGC
    plot <- ggplot( data = dat_plot2[dat_plot2$Type=="IGC vs DGC",] , aes( x = Class_num , y = MutRate , fill = Class ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Mutation Rate")+
    theme(panel.grid = element_blank())+
    scale_fill_manual(values=col) +
    ylim(0,1.05)+
    geom_text(aes(label=p_text , y = 1 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=4 , face='bold' , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='none',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 8,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.ticks.x = element_blank(),
                axis.text.x = element_text(size = 8,color="black",face='bold') ,
                axis.line = element_line(size = 0.5))

    out_name <- paste0( images_path , "/MutRate_" , gene , ".IGC_DGC.pdf" )
    ggsave( out_name , plot , width = 2.5 , height = 5 )


    ## IM vs GC
    plot <- ggplot( data = dat_plot2[dat_plot2$Type=="IM vs GC",] , aes( x = Class_num , y = MutRate , fill = Class ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Mutation Rate")+
    theme(panel.grid = element_blank())+
    scale_fill_manual(values=col) +
    ylim(0,1.05)+
    geom_text(aes(label=p_text , y = 1 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=3 , face='bold' , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='none',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 12,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.ticks.x = element_blank(),
                axis.text.x = element_text(size = 8,color="black",face='bold') ,
                axis.line = element_line(size = 0.5))

    out_name <- paste0( images_path , "/MutRate_" , gene , ".IM_GC.pdf" )
    ggsave( out_name , plot , width = 2.5 , height = 5 )
}

## 不同来源
for( gene in unique(dat_plot$Hugo_Symbol) ){
    dat_plot2 <- subset( dat_plot , Hugo_Symbol == gene & From != "All" )
    dat_plot2$Type <- ""
    dat_plot2$Type[dat_plot2$Class %in% c("IM" , "GC")] <- "IM vs GC"
    dat_plot2$Type[dat_plot2$Class %in% c("IGC" , "DGC")] <- "IGC vs DGC"
    dat_plot2$Type <- factor( dat_plot2$Type , levels = c("IM vs GC" , "IGC vs DGC") )
    dat_plot2 <- subset( dat_plot2 , From != "Utokyo" )
    dat_plot2$From <- factor( dat_plot2$From , levels = c("NJMU" , "OncoSG" , "TCGA" , "TMUCIH" , "Utokyo") )  
    dat_plot2$Class_num <- paste0( dat_plot2$Class , "(" , dat_plot2$SampleNum , ")" )  
    dat_plot2$Class_num <- factor( dat_plot2$Class_num , levels = unique(dat_plot2$Class_num)[order(unique(dat_plot2$Class_num) , decreasing=T)] )  
    
    ###########################################################################################
    ## IGC vs DGC
    plot <- ggplot( data = dat_plot2[dat_plot2$Type=="IGC vs DGC",] , aes( x = Class_num , y = MutRate , fill = Class ))+
    geom_bar(position = "stack", stat = "identity") + 
    theme_bw()+
    labs(x="",y="Mutation Rate")+
    facet_grid(.~From,scales="free")+
    theme(panel.grid = element_blank())+
    scale_fill_manual(values=col) +
    ylim(0,1.05)+
    geom_text(aes(label=p_text , y = 1 ,x = 1.5),parse = TRUE,size=4)+
    geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=4 , face='bold' , color="black")+
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
                legend.position ='none',
                legend.title = element_blank() ,
                panel.grid.major=element_line(colour=NA),
                plot.title = element_text(size = 12,color="black",face='bold'),
                legend.text = element_text(size = 12,color="black",face='bold'),
                axis.text.y = element_text(size = 8,color="black",face='bold'),
                axis.title.x = element_text(size = 12,color="black",face='bold'),
                axis.title.y = element_text(size = 12,color="black",face='bold'),
                strip.text.x = element_text(size = 7 , face = 'bold'),
                axis.ticks.x = element_blank(),
                axis.text.x = element_text(size = 8,color="black",face='bold') ,
                axis.line = element_line(size = 0.5))

    out_name <- paste0( images_path , "/MutRate_" , gene , ".IGC_DGC.From.pdf" )
    ggsave( out_name , plot , width = 5 , height = 5 )
}